Spectrum of low energy excitations in the vortex state: 
comparison of Doppler shift method to quasiclassical approach 
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We present a detailed comparison of numerical solutions of the quasiclassical Eilenberger equations 
with several approximation schemes for the density of states of s- and d-wave superconductors in 
the vortex state, which have been used recently. In particular, we critically examine the use of 
the Doppler shift method, which has been claimed to give good results for d-wave superconductors. 
Studying the single vortex case we show that there are important contributions coming from core 
states, which extend far from the vortex cores into the nodal directions and are not present in the 
Doppler shift method, but significantly affect the density of states at low energies. This leads to 
sizeable corrections to Volovik's law, which we expect to be sensitive to impurity scattering. For 
a vortex lattice we also show comparisons with the method due to Brandt, Pesch, and Tewordt 
and an approximate analytical method, generalizing a method due to Pesch. These are high field 
approximations strictly valid close to the upper critical field B C 2- At low energies the approximate 
analytical method turns out to give impressively good results over a broad field range and we 
recommend the use of this method for studies of the vortex state at not too low magnetic fields. 

PACS numbers: 74.20.Rp, 74.60.-w, 74.25.Bt 



I. INTRODUCTION 



There is now a growing number of candidate systems 
for unconventional or strongly anisotropic superconduc- 
tivity. Besides the high-T c cuprates, which are believed 
to be d-wave superconductors, there are indications f< 
unconventional superconductivity in Sr2Ru04 J!-.UPt3 
an organic superconductor ,0 the «-(ET)a saltsn and a 
possible unconventional, ferromagnetically driven super- 
conductivity is discussed in the recently found supercon- 
ducting ferromagnets ZrZn2, URhGe and UGe2.^-Also, 
in MgB-2 phonon,-piediated, strongly anisotropicoLl or 
two-gap scenariosel have been proposed. In such sys- 
tems the strong momentum dependence of the super- 
conducting order parameter can lead to interesting new 
behavior, especially if there are gap nodes present at 
the Fermi surface. In particular, the low-energy ex- 
citations in the vortex state of such systems are ex- 
pected to display.. unronyrnlional hohavipr, as has been 
studied recently.MBBE^BBBBB In these stud- 
ies the so-j«aliedirO©Bpler shift method has been used 
frequently.BBBBEJ In this method, which dates back 
to the early day*, of the theoretical study of type II 
superconductors,E2l the quasiparticle excitation energy is 
approximated by taking into account only the Doppler 
shift due to the local supercurrent flow. While this 
method neglects vortex core scattering and certainly is 
a bad approximation for conventional s-wave supercon- 
ductors, Volovik has argued that in the case of super- 
conductors having gap nodes, the low energy excitations 
are dominated by the extended quasiparticle states out- 
side the vortex core and the Doppler shift method thus 
should give much better results, at least for high-K super* 
conductors at sufficiently low external magnetic fields .lI 
Based on the Doppler shift method, Volovik predicted 



that the density of states at the Fermi level in the vortex 
state of a superconductor with line nodes should vary as 
the square root of the magnetic field, instead of the lin- 
ear variation expected for a conventional superconductor. 
Indeed, low-temperature specific heat measurements on 
high-T c cuprates have bp pp sbrpnp to be consistent with 
such a field dependenceOBBBB 

However, the Doppler shift method is not a rigorous 
quasiclassical approximation in the sense of an expan- 
sion in terms of (/cf£) _1 , where 1zf is the Fermi momen- 
tum and £ the coherence length. Instead, the appropri- 
ate method is the solution of tha-Fil<?riberger equations 
for the quasiclassical propagator .BBB Within this ap- 
proach, the contribution of vortex core states and vortex 
core scattering is included. Of course, the Eilenberger 
approach requires the solution of a set of transport equa- 
tions, which is more involved than the Doppler shift cal- 
culations. 

In the present work we want to present detailed com- 
parisons of Doppler shift calculations with solutions of 
the Eilenberger equations for s- and d-wave superconduc- 
tors, in order to clarify where the Doppler shift method 
is a good approximation and where it fails to give quan- 
titative agreement with the fully quasiclassical solution 
of the Eilenberger equations. We are also going to dis- 
cuss two other approximate methods, which are expected 
to be good at magnetic fields close to the upper critical 
field B C 2- Indeed, in this field range it turns out that an- 
alytical results for the density of states averaged over a 
unit cell of the vortex lattice can be obtained, which are 
superior in accuracy than the Doppler shift results. A nu- 
merical solution of the Eilenberger equations is simplified 
considerably due to a recently found parametrization, 
which transforms the Eilenberger equations into scalar 
differential equations of the Riccati type. 
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In the next section we want to briefly present the Eilcn- 
berger approach and its mapping to scalar Riccati equa- 
tions. We will outline, how the Doppler shift method can 
be obtained from the Riccati equations. Section 



III 



is de- 
voted to the study of the local density of states for the 
single vortex, while in section IV we are going to present 
comparisons of the density of states for the vortex lattice. 



II. EILENBERGER APPROACH 

For a layered system with cylindrical Fermi surface 
and a spin-singlet superconducting order parameter the 
Eilenberger equation for the normal and anomalous 
com pjQjjjen ts g and / of the quasiclassical propagator 



(ie n H — vf - A) +/7/r'/. • V 



f(f, 0, ie n ) = 
2«Kr,e,z6 n )A (r,e) (1) 



Here, Vp is the Fermi velocity pointing into the direction 
0, A is the vector potential due to the internal magnetic 
field within the system, e n are the fermionic Matsubara 
frequencies, and A (r, 0) denotes the spatially varying 
and momentum dependent order parameter. Eq. (Q) has 
to be supplemented by a normalization condition, which 
in our case readscJ 

{g(r,Q,ie n )} 2 + f(r,Q 1 ie n )f*(r,e + ir,ie n ) = 1 (2) 

The normalized local density of states N(r, e) is obtained 
from g after an analytic continuation to real frequencies 
and an angular average over the Fermi surface: 



N(f,e) 



1 

27 



dO Re {g(r,B, 



ic. 



It has been shown in Refs. [h],^8| that Eqs. ([!]) and (|^) 

can be mapped onto scalar Riccati equations along real 
space trajectories r(x) running parallel to the direction of 
the Fermi velocity, if one introduces two scalar complex 
quantities a(x) and b(x) 



l(f(x)) 



1 — a(x)b{x) 



f(r(x)) 



2ia(x) 



(4) 



l + a(x)b(x) l + a(x)b(x) 

where a(x) and b(x) obey the following Riccati equations 
d 

hv F —a(x) + \2e n (x) + AUx)a(x)] a(x) - A(x) = 
ox L J 

(5) 

hv F -^-b(x)~ \2eJx) + A(x)b(x)]b(x) + AUx) = 
ox 

(6) 

Here, ie n (x) = ie n + -vp ■ A{x). The initial values for the 
two quantities a(x) and b(x) in the bulk superconductor 



have to be taken as 



a(— oo) = 



A(-oo) 



&(+oo) 



:„ + + |A(-oo)| 2 

At(+oo) 
, l + v /^ + |A(+^)| 2 



(7) 



(8) 



Once the order parameter field Ao(r, 0) and the vec- 
tor potential A(r) are known, the initial value prob- 
lem for the scalar differential equations (H) and ((6) can 
be solved by standard numerical techniques (adaptive 
Runge-Kutta method for example). For a homogeneous 
bulk superconductor one confirms that Eqs. (|t]) and (||) 
indeed fulfil Eqs. (|) and (^). In order to find the local 
density of states Eq. (j^) for a given point r and energy e 
it is necessary to solve the Riccati equations (||) and (||) 
for a bundle of trajectories running through the point r 
with different angles 0. 

It is instructive to see how the Doppler shift method 
can be obtained from Eq. (^|). For that purpose we first 
decompose the order parameter A (a;) into amplitude and 
phase via 

A(x) = A(x)e l ^ x) (9) 
Introducing the function 

a(x) = a^e-^ 
one arrives at an equation for a: 
d 



(10) 



Hvf-^-o,(x) 
Ox 



[2e„ + 2ivp ■ mv s (x) + 
+A(x)a{x)] a(x) - A(x) = (11) 



where v s {x) — ^ yfiV^(x) — ^A(x)j is nothing but 

the gauge invariant superfluid velocity of the supercur- 
rent field distribution. If we now neglect the spatial 



(12) 



i0 + )} (3) derivative -^a(x) we find from Eq. ( p.1 



a{x) 



A(x) 



e n {x) + ^e n {x) 2 + A(xf 



with e n {x) = e n + imvp ■ v s (x). Thus, we rediscover the 
bulk result apart from a Doppler shift in energy. If we set 
A (a;) equal to its homogeneous bulk value using Eqs. (H), 
and (0) we just obtain the usual Doppler shift equationEl 



(13) 



2jv 



— d0Re 



1 

2~ix 



|e — mvp ■ v s (r)\ 



^ {e - rnvp ■ v s {r)f - |A(0)| 2 



From this derivation we learn that the Doppler shift 
method neglects the gradient term in the Eilenberger 
equations. This neglection is expected to be a reason- 
able approximation as long as the Doppler shift energy 
mvp ■ v s (r) is small compared to the local gap energy 
A(r, 0). As is well known, this approximation fails close 
to the vortex cores, where the superfluid velocity di- 
verges. However, as we will see later, this approximation 
also fails in the vicinity of gap nodes of A(r, 0). 
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III. THE SINGLE PHASE VORTEX 

In order to make quantitative comparisons between the 
Doppler shift method and the Eilenberger approach, we 
first want to study the single vortex case. In particular, 
we restrict ourselves to the pure 'phase vortex' for which 
the amplitude of the order parameter A(0) is assumed to 
be constant as a function of r and only its phase is vary- 
ing. For a d-wave superconductor taking the magnetic 
field along the c-axis direction we then have 



A(f, 9) = A cos(26)e 



i<f> 



(14) 



where <fi denotes the polar angle of r. We choose this 
model for the vortex, because it should be the 'best case' 
for the Doppler shift method, since it neglects any spa- 
tial variation of the amplitude of the gap. Such a model 
for the vortex is expected to be reasonable for a high- 
re superconductor at low magnetic fields of the order of 
the lower critical field H c i. Although this phase vortex 
does not possess a core in the usual sense, we will see 
that there still exists a core region in the Eilenberger 
approach, which is due to the gradient term. One ob- 
serves that the gradient term introduces a length scale 
£ = hvp/A into the problem, even if the amplitude A 
itself does not vary. For the single vortex the Doppler 
shift energy diverges like r -1 and we have 



mv p ■ v s (r ) 



2r 



~27 



sin(6 - (f) (15) 



where is the unit vector in (^-direction. This equation 
holds as long as r is smaller than the penetration depth. 
Otherwise the screening of the magnetic field has to be 
taken into account as well. 

We have calculated the local density of states for the 
Doppler shift method using Eqs. (|l|) and (|f). The 
local density of states for the Eilenberger approach is 
obtained from Eqs. ([|) and (||) numerically using the 
Riccati method outlined above. In order to facilitate the 
solution a small imaginary part 5 < O.OlAo has been 
added to the energies on the real axis. 

In Figs. and E we are showing our results for the 
energy dependence of the local density of states for a d- 
wave superconducting state. Here and in the following 
the density of states is normalized to the normal state 
value. In Fig. 0(a) and (b) the results in an angular 
direction of <p = are shown for distances r = 1 and 
r = 3 from the vortex center (in units of the coherence 
length £). Figs. ||(a) and (b) show the corresponding 
results for an angular direction of <fi — 7r/4 (nodal direc- 
tion). As becomes clear from these plots, the Doppler 
shift method gives resonable results at sufficient high en- 
ergy and also becomes better at distances further away 
from the vortex center. However, at distances from the 
vortex of the order of the coherence length the Doppler 
shift fails to give the position of the peaks correctly, 
which result .from core state contributions or scattering 
resonances .E3E3 The satellite peaks in the Doppler-shift 
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FIG. 1: Local density of states in the vicinity of a d-wave 
phase vortex at distances (a) r = 1 and (b) r — 3 (in units of 
the coherence length) from the vortex center in the antinodal 
direction = 0. The dashed line shows the result from the 
Doppler shift calculation, while the solid line is the result from 
a numerical solution of the Eilenberger equations. 



method seen in Figs. |l| and appear at the Doppler 
shifted gap energy, as has been discussed recently in Ref. 
|l5| . A discussion of the position of the peaks for the d- 
wave vortex within the quasiclassical approximation can 
be found in Refs. |l(i| , p9| . 

It has been noted earlier that in. a c?-wave vortex 
the core states are not truly localizedElJH Nevertheless, 
they give important contributions to the local density 
of states, which are not captured by the Doppler shift 
method, but are present in the numerical solution of 
Eilcnbcrger's equations. We emphasize that these res- 
onances are still present, although our phase vortex does 
not possess a variation of the magnitude of the gap. 
Thus, they cannot easily be understood as 'bound states' 
in the potential well of the gap at the vortex core. In- 
stead, these resonances arise due to scattering of the 
quasiparticles at the phase gradient around the vortex 
center and thus can be interpreted as Andreev scatter- 
ing resonances. For comparison, we also did calculations 
where we multiplied Eq. (|l4|) by tanh(r/£ c ), where £ c is 



4 




FIG. 3: Zero energy local density of states at a distance of 
r = 5 from the vortex center plotted as a function of polar 
angle <f>. The circles are numerical results obtained from the 
solution of Eilenberger's equations using different imaginary 
parts: S — 0.001 (full circles) and S — 0.01 (open circles). 
The dashed line is the result obtained from the Doppler shift 
(DS) method. The solid line shows the angular dependence 
expected from a solution of the Bogoliubov-de Gennes (BdG4 
equations linearized around the gap nodes due to Mel'nikovEj 




o 1 1 

0.5 1 1.5 2 

6 A, 

FIG. 2: Same as Fig. |l], but in the nodal direction <f> = n/4 
from the vortex center. 



an adjustable parameter independent of £, which allows 
us to smoothly cross over from a vortex possessing a core 
in the usual sense to our phase vortex, taking the limit 
£ c — > 0. These calculations confirm that the core state 
contributions do not disappear in this limit, but instead 
smoothly evolve into the peaks seen in Figs. | and §. 

Even some distance away from the vortex center, the 
Doppler shift does not give the correct behavior at low 
energies, however. As can be seen in Fig. |(a) and (b), 
there is a small peak at low energies in the local density of 
states in the quasiclassical solution. This peak is mainly 
visible in the vicinity ofJJae nodal direction <fr = 7r/4, 
as has been noted beforepj and is due to extended core 
states, which can 'leak out' of the vortex core region due 
to the nodes in the gap function. This effect is also not 
present in the Doppler shift calculation. 

This becomes more transparent from Fig. ^, where we 
show the zero energy local density of states at a distance 
r = 5 from the vortex center as a function of polar an- 
gle 4>. The symbols show the results from our numerical 
solution of Eilenberger's equations, open circles for an 
imaginary part of 8 = 0.01 and solid circles for 5 = 0.001 
(in units of Ao). The dashed line shows the Doppler 
shift result. (A discussion of the solid line is found be- 



low). From this plot it becomes clear that the solution 
of Eilenberger's equations yields increased contributions 
to the local density of states at zero energy especially 
close to the nodal direction <\> = 7r/4, which the Doppler 
shift method is not able to capture. We can trace back 
these increased contributions to quasiparticle trajecto- 
ries passing near by the vortex center with a momen- 
tum direction close to the nodal direction O = 7r/4 by 
looking at the momentum resolved local density of states 
N(r, <p, 0, e) = Re {g(r, <p, 0, e)}. On these particular 
trajectories the gap is small and the corresponding wave- 
function of the quasiparticle extends very far, 'leaking 
out' of the core region. These nonlocal effects are missed 
by the Doppler shift method, since it only knows about 
the local supercurrent flow. As Fig. || shows, this effect 
strongly depends on the imaginary part S chosen. It is 
clear that 6 introduces a finite mean free path into the 
problem, which limits the extend to which these quasi- 
particle states can contribute to the local density of states 
far from the vortex core. Thus, one should expect that 
this effect is sensitive to impurity scattering. 

One notes that in Fig. |^ the results from Eilenberger's 
equations display peaks slightly off the <f> = 7r/4 direction. 
The position of the peaks depends on 5 and moves closer 
to 7r/4 when 5 is reduced. The suppression directly at 
7r/4 is due to the presence of the d-wave gap node in 
this direction, because a trajectory passing through the 
vortex center in this direction does not 'see' a gap and 
thus no resonant Andreev scattering takes place. (For 
comparison see also Fig. 8d in Ref. ||^). We observe 
that the main contribution to the local density of states 



5 



o 





— 6 


i i i i 1 i i i i 1 i i ,j 

= o.oi : 




— <5 


= 0.005 , 


'- 


— 6 


= 0.0025 : 










0.05 0.1 0.15 0.2 

1/R 

FIG. 4: Average density of states within a circle of radius R 
(in units of the coherence length £) around the vortex plotted 
as a function of 1/7?. The dashed line is Volovik's result 
obtained from the Doppler shift method. The symbols are 
numerical results obtained from a solution of Eilenberger's 
equations using different imaginary parts: 8 = 0.01 (open 
circles), 8 = 0.005 (full circles), and 8 = 0.0025 (full squares). 
Here, 8 is measured in units of the gap amplitude Aq. 



in this case is coming from trajectories slightly off the 7t/4 
direction, depending on the momentum width induced by 
5. 

This effect also has important consequences for the 
magnetic field dependence of the density of states at zero 
energy. In Fig. [I| we show the density of states aver- 
aged over a circle of radius R around the vortex center 
as a function of 1/R, which at low fields and for high 
k superconductors is proportional to the square-root of 
the magnetic field. The Doppler shift method yields a 
square-root dependence of the density of states as a func- 
tion of magnetic field (a linear variation as a function of 
1/R), which has been noted first by VolovikB However, 
the quasiclassical solution leads to important deviations 
from this law. First of all, the slope of these curves is 
much higher, and secondly this effect is very sensitive 
to the imaginary part chosen in the calculation, as is 
clear from the discussion above. While we believe that 
the field dependence of these curves will probably not be 
distinguishable experimentally from the square-root field 
dependence, a systematic study of the influence of impu- 
rity content on the field dependence of the specific heat 
at low temperatures might be able to detect the sensi- 
tivity of the slope to impurity scattering and would be a 
valuable confirmation of the presence of these extended 
core states. 

Since we are working here with a quasiclassical ap- 
proximation one might ask, how important quantum ef- 
fects like the Aharanov-Bohm effect might be on these 
extended core states. Indeed, the Aharanov-Bohm ef- 
fect for a (i-wave vortex has been studied recently by 
Mel'nikov within a solution of the Bogoliubav-de Gennes 
equations linearized around the gap nodes £3 Within this 



quantum mechanical calculation Mel'nikov showed that 
the angular dependence of the zero energy density of 
states far from the vortex core is determined by the Dirac 
cone anisotropy a = vp/va, where va is the quasi- 
particle velocity tangential to the Fermi surface at the 
node, which is given by the slope of the gap (in our case 
Va = 2Ao/hkp and thus a = kp^/2). For_the angular 
dependence Mel'nikov found the expressions 

N(r,4>- ~,e = o) = (16) 



n r I \J a 2 cos 2 4> + sin 2 4> y a 2 sin 2 cj> + cos 2 4> \ 

(Note, that in Ref. [32] the nodal direction corresponds 
to tj> = 0. Therefore, we are shifting this expression 
here by 7r/4. The prefactor A/ir can be found from 
Mel'nikov's result normalizing it to the normal density 
of states). For comparison this angular dependence is 
shown in Fig. || for a = 20 as the solid line. The angu- 
lar dependence becomes more pronounced for higher val- 
ues of a, but estimates from thermal conductivity mea- 
surements on high-T c euprates indicate anisotropies in 
the range a ~ 15 — 20. Clearly, the contribution from 
the extended core states is contained in both the quan- 
tum mechanical and the quasiclassical calculation. Note, 
that the full quantum mechanical information about the 
particle-hole coherence along the classical trajectories is 
contained in the quasiclassical approach (while the quan- 
tization perpendicular to the trajectories is neglected). 
Our quasiclassical approach corresponds to taking the 
limit 2a = kp£ — * oo, but does not depend on a lin- 
earization around the gap node. From the comparison 
in Fig. [| we expect that full quantum interference, i.e. 
interference of quasiparticles running classically different 
ways around the vortex line, will limit (via the Dirac 
cone anisotropy) the angular dependence of the density 
of states. But this comparison also shows that a very pure 
sample is needed, if one wants to observe these quantum 
mechanical effects, since at shorter mean free paths the 
angular dependence is expected to be limited by the im- 
purity scattering rate. 

IV. VORTEX LATTICE 

We now want to focus our attention on the density of 
states in the vortex lattice. In the vortex lattice, the or- 
der parameter profile as well as the superfluid flow around 
the vortices has to be modified as compared to the sin- 
gle vortex case in Eqs. ( fTi]) and ([To]). Specifically, here 
we will use Abrikosov's wavefunction, which is the exact 
solution close to the upper critical field B C 2- An arbi- 
trary vortex lattice A can be spanned by two lattice vec- 
tors Ri — (wi,0) and R2 = (Re W2,Im W2), where we 
combined the x- and y-component of R2 into the com- 
plex quantity C02 and we have chosen the x-axis into the 
R\ direction. Without loss of generality we may assume 
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u>i > and Im u>2 > 0. Using the quantities oji and cl>2 we 
can generalize Abrikosov's wavefunction in the following 
form: 



-0a (x,y) 



1 



\ " 



exp 



7T (ixy - y 2 ) 
uji Im 



(2n + 1) . w 2 , 1N 
(x + ly) + itt — n (n + I) 



(17) 



where M = 



Im uJ2 

0)1 



1/4 



is a normalization 



2im^ ex P 

factor such that \4>a\ 2 averaged over a unit cell Ca of the 
vortex lattice becomes unity. For a square lattice we have 
LU2 = while for a triangular lattice ui2 — 1+ 2 ^l- 
The average magnetic field B av in the superconductor 
points into the z-direction and the gauge has been chosen 
such that the vector potential obeys the relation 



.4 



1 



-fX Br, 



(18) 



Due to the flux quantization the average magnetic field 
is related to the area of the unit cell of the vortex lattice 
yielding the relation B av = 2 —- h £ - . The superfluid 
velocity field v s (r) for the Doppler shift calculation can 
then be obtained from Eq. ([17|) using the relation 



v s (r) = 



Ami 



Im 



-.4 



(19) 



Due to the n 2 term in the exponent the sum in Eq. (|]J 
quickly converges and only a few terms are sufficient to 
find ip\ [x,y). For our numerical solution of the Eilen- 
berger equations for a d-wave superconductor we then 
use the order parameter 



A(f,6) = A cos(29)VA(fO 



(20) 



When solving the Riccati equations Eqs. (j^) and (^|) for 
the vortex lattice we have to ensure that a solution peri- 
odic in the lattice is found. We achieve that by utilizing 
our small imaginary part 8 for the energies. For a given 
point in space and momentum we choose the starting 
point of the corresponding trajectory a long distance of 
order 1/8 unit cells away. This gives the solution a dis- 
tance to relax to the periodic solution sought, where the 
relaxation is provided by 8. We have explicitly checked 
that the result is periodic and doesn't change anymore, 
if we repeat the calculation with a longer trajectory. 



A. An approximate analytical solution 

At high magnetic fields close to the upper critical 
field B C 2 it is possible to obtain analytic results for 



the density of states averaged over a unit cell of the 
vortex lattice, generalizing a method due to Pesch,E£l 
which has been used cecently also for the study of d- 
wave superconductors .E3 These analytic results can be 
obtained as follows. The first step is to replace the func- 
tion g on the right hand side in Eq. ([!]) by its value aver- 
aged over a unit cell of the vortex lattice. This approxi- 
mation is certainly valid close to the upper critical field, 
because there g does not vary strongly within a unit cell. 
Then we have 

Lf(r, 9, ien) = 2(g(f, 6, ie n ))c A A(r, 9) (21) 

where (• • • )c A denotes an average over a unit cell of r 
and the operator L is given by 

L = L (6) = 2 (e n - i-v F (6) • A) + M F (6) • V (22) 

Assuming without loss of generality that eri-n* Eq. (|2l| ) 
can be inverted using the operator identitycj 



L~ x = 



ds 



(23) 



Since {g)cx does not depend on the variable r anymore, 
L only acts upon A and we find 



f(f, 6, ie n ) = 2(g(f, 9, ie n ))c A 



ds e- sL A(f,9) 

(24) 

. (H) we can calcu- 



Using the normalization condition Eq. 
late the cell average of the square of g: 

(g 2 (r,<d,ie n ))c A 
= 1 - (/(r, 9, ie n )f*(r, 9 + tt, ie n ))c A 
= l-(g(r,<3,ie n )) 2 CA P A (<d,ie n ) (25) 



where 
PA{Q,ie n ) 



ds- 



dsj 



ii Jo 



((e- s + L A(f,9)) (e- 



(26) 

At(f,9 + 7r))) OA 

Here we have used the relation <^(9 + tt) = g(9) and 
inversion symmetry L^(Q + tt) — L(Q). Close to the 
upper critical field we may assume as a second step that 



(g 2 (r,Q,ie n ))c A = {g{?,<d,ie n ))c A 



(27) 



which amounts to neglecting spacial fluctuations of g. 
Using this approximation Eq. (25 ) becomes a closed equa- 



tion for (g)c A and we finally find 

(g(r,Q,ie n ))c A = 



1 



v /l + P A (9,ie„) 



(28) 



As it turns out, Eq. (|26|) can be evaluated analytically 
for the_Abrikosov vortex lattice Eq. (|i~7|). After some 
algebraE3 we find 



Pa(9, 1£ „) 



2 A 2 cos 2 (29) 



z 2 [l - w {iz)\ (29) 
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where 



(30) 



Here, the function w is Dawson's integral and is related 
to the complement of the Error function via 



w{iz) 



1 

in 



'dt 



t 



erfc(z) 



(31) 



Some properties of this function can be found in Ref. 
|37| . We can make these equations a little bit more trans- 
parent, if we introduce two length scales: the coherence 
length £ = vph/A Q and a\ = Wilm u> 2 = &o/B av with 
$0 being the flux quantum. a\ is the area of a unit cell 
of the vortex lattice and for a square lattice oa is just the 
distance between neighboring vortices. We remark that 
also the coherence length £ depends on magnetic field, 
since Ao has to be determined from the gap equation in 
the presence of the field. Expressed in these quantities 
we have 



2 Qa £» 
7T £ A 



(32) 



To obtain the cell average of the density of states, 
Eq. ( p8| ) has to be integrated over the angle 6 and an- 
alytically continued to real frequencies ie n — > e + i0 + . 
Since Pa only depends on 9 via the cos 2 (26) term, the 
angular integral is just a complete elliptical integral and 
we find the analytical result 



1 f 2n 

N(e) = — J d6Re( ff (f,e, 



+ i0 + )) 



C A 



-Re{K{k)} 

7T 



with 



,2__4a A 

k - 




(33) 



(34) 



For comparison, for an s-wave superconductor the 
cos 2 (20) factor in Eq. (|2^ ) has to be dropped and we 
just find 



N(e) = Re 



1 



VI - fc 2 



(35) 



It is instructive to consider the low and high field lim- 
its of these expressions. At high magnetic field B C 2 the 
coherence length £ diverges. In this limit k — > and we 
find the normal state result N(e) — ► 1 for both s- and 
e?-wave superconductor as one should expect. For low 
magnetic fields £ becomes constant, but oa diverges. In 
this limit onC|-C|an use the asymptotic expansion of Daw- 
son's integraltll 



_ . . 1 

1 - ^TTZ W (lZ) ~ 



(36) 




FIG. 5: Average density of states in the vortex state of a 
d-wave superconductor calculated using the analytical for- 
mula Eq. ([33]) for different field strengths: S av = (solid 
line), Bav = 0.1B C 2 (dashed line), B av — 0.3-B C 2 (dashed- 
dotted line), and _B av = 0.5B C 2 (dotted line). The energy 
is normalized to the gap value at zero field and temperature 
Aoo = A (T = 0,B = 0). The inset shows the field depen- 
dence of the zero temperature gap Ao/Aqo- 



Then we have k 2 ~ Ap/e 2 . Inserting this into Eqs. (|32^ 
and (^5|) we just rediscover the zero field density of states 
of the d- and s-wave superconductor, respectively. It is 
interesting to note that this high field approximation also 
appears to give the correct zero field limit. 

In Fig. H we show the density of states calculated from 
Eq. (|3^) for different magnetic fields at zero temperature. 
In order to find the appropriate value of <2a/£ for a given 
field and temperature, it is necessary to solve the gap 
equation in order to find Ao(T, B av ). Using the approxi- 
mation Eq. (^4]) the gap equation for the <i-wave case can 
be cast into the form 



In 



T 
% 



(37) 



^200^(20) 



I~~kz w (iz) 



v /l + P A (e,ie„) 



with z depending on e n and B av via Eq. ([50|). The an- 
gular integration over O can even be further reduced to 
complete elliptical integrals, if desired. The field depen- 
dence of the gap at zero temperature obtained from this 
equation is shown in the inset to Fig. |^. The energy scale 
in Fig. |^ has been normalized to the zero field and tem- 
perature value of the gap A o = A (T = 0, B av — 0). In 
the inset it can be noted that the gap Ao as a function of 
magnetic field rises above Aoo at low field before reach- 
ing its zero field value Aoo at B av = 0. This is a property 
of Eq. (37) and we attribute it to the fact that the ap- 
proximation used in the present section is valid only at 
high field and at zero field. Fig. [| shows that, although 
the gap closes as a function of field B av , the peaks in the 
density of states move upwards in energy. This means 
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that in the vortex state the peak-to-peak distance in the 
density of states is not an appropriate measure of the en- 
ergy gap Ao anymore. This fact is important to realize 
in interpretation of tunneling spectra in the vortex state, 
for example. 



B. The approximation due to Brandt, Pesch, and 
Tewordt 

In the literature there exists another approximation for 
the average density of states in type II superconductors 
at high magnetic fields, which is due to Brandt, Pesch, 
and TewordtEa and also has beeia-j*sed recently to study 
unconventional superconductors J23e3 In their original pa- 
per these— authors derived their method from Gorkov's 
equation.^ We discovered a method to re-derive their 
results from the Eilenberger equation using the Riccati 
equations (|^) and (Q), which we want to briefly sketch 
here. 

Instead of replacing the function g by its average value 
in Eq. (|l|) we can use an approximation similar in spirit 
in the Riccati equations (||) and (||) , linearizing them by 
replacing the terms A^a and Ab by their averages over a 
unit cell. For the function a one obtains the equation 



La(r,Q,ie n )=A(f,G) 



(38) 



where L is just the operator in Eq. (|22|), except that e„ 
now has to be replaced by 



£« = e n + ^(A^a)c A 



(39) 



Eq. (|3S|) can be inverted the same way as above and after 
some algebra we obtain for the average (A'a)c A - 



<At a ) CA = / ds (At(r, 6)e- s£ A(f, 0)>c A 



f'A 



= Ao cos (2Q)—^w(iz) — 2 (e„ — e„ 



(40) 



where now 



y/nvph ~ V 7T £ A 

Eq. ( ^o|) is an implicit equation for e„ or equivalently for 
z and can be written in the form 

£^+W(20)-|L»(«) (42) 

From the inversion of Eq. ( |3^ ) and the equivalent one for 
the function b one can also obtain the cell average of the 
product ab: 

(ab) Ca = cos 2 (26) [1 - V^z w (iz)] (43) 



Using this result we can obtain an approximation for the 
cell average of the propagator g using Eq. M) 



(g(f,Q,ie n ))c A 



1 + {ab)c A 
2 



1 



- 1 (44) 



1 + cos 2 (20) [1 ~ V^z w (iz)] 

where z has to be obtained from a solution of Eq. (42) 
for each e„ and 0. The set of equations ( |42] ) and (44) 
just corresponds to Eqs. (18) and (16) in Ref. Ejl. (The 
parameter A in that work corresponds to a\/V2ir in our 
case). 

This derivation shows that the approximation due to 
Brandt, Pesch, and Tewordt can also be justified from 
quasiclassical Eilenberger theory. However, in contrast to 
the approximate analytical method derived in the previ- 
ous subsection, here one has to solve a nonlinear implicit 
equation, which makes the calculation of the density of 
states more difficult. An expansion of Eqs. (44) and ( psj ) 
in terms of Ao shows that both methods give the same 
result up to order |A | 3 . 



C. Comparison with numerical results 

We now want to present a comparison of the four meth- 
ods outlined above: the Doppler shift method (DS), the 
approximate analytical solution ( AA) , the method due to 
Brandt, Pesch, and Tewordt (BPT), and a full numerical 
solution of Eilenberger's equations (EE) using the Riccati 
equations. In all four cases we will base our calculations 
on an Abrikosov square lattice. However, the results do 
not differ very much, if a triangular lattice is used. In 
fact, the density of states within the methods AA and 
BPT does not depend on the lattice structure, as is clear 
from the derivation above. 

In the following we will compare our results using the 
parameter £/a\. This parameter is a measure of the av- 
erage magnetic field B av inside the superconductor and 
at low field strengths £,/a\ oc s/B av . At higher field 
strengths its field dependence has to be determined from 
a solution of the gap equation, because £ depends on Ao, 
which decreases with increasing field. From such a so- 
lution we have determined that £/a A = 1 corresponds 
to about -B av ~ 0.5B C 2 and £/a\ = 0.2 corresponds to 
about B av w 0.04i? C 2 at zero temperature. 

In Fig. H we show the cell average of the density of 
states as a function of energy for the four methods. In 
Fig. ||(a) the comparison for £/a\ = 1 is shown, Fig. ||(b) 
shows the results for £/a\ = 0.2. At high magnetic 
field (£/<xa = 1) the numerical solution of the Eilen- 
berger equations (EE, black dots), and the methods AA 
(solid line) and BPT (dotted line) are all in very close 
agreement with each other, the approximate analytical 
(AA) result being somewhat closer to the numerical so- 
lution. The Doppler shift method (dashed line), how- 
ever, strongly deviates from these results. This result is 
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FIG. 6: Energy dependence of the density of states averaged 
over a unit cell of the vortex lattice for (a) £,/cla = 1 and (b) 
£/oa = 0.2, corresponding to magnetic fields of about 0.5B C 2 
and 0.04B C 2, respectively. The dots are numerical solutions of 
the Eilenberger equations, the dashed line shows the Doppler 
shift result, the dotted line the result from the method due 
to Brandt et al, and the solid line the approximate analytical 
result described in the text. 



FIG. 7: Zero energy density of states as a function of the 
parameter £/aA for an (a) d-wave and an (b) s-wave super- 
conductor. The four methods are denoted by the same line 
patterns as in Fig. ^. 



not surprising, because the Doppler shift method neglects 
any contributions from vortex core scattering, while the 
methods AA and BPT are expansions around B C 2 and 
are expected to give better results at high fields. 

Fig. ^(b) shows the comparison at a considerably 
smaller field (£/a\ = 0.2). In this field range the four 
methods yield different results, especially close to the 
gap edge. Xhe numerical result clearly shows Tomasch 
resonances,c3 which the other three methods are not 
able to get. At high energy all methods converge to 
each other. At low energies, which are especially impor- 
tant for the thermodynamics of the system, the Doppler 
shift method gives the strongest deviation, while the AA 
method is closest to the numerical solution. 

In Fig. we show the density of states at zero energy 
as a function of the parameter £,/a\ for both d- and s- 
wave superconductor. From these plots it can be seen 
that in this zero energy limit the AA method appears to 



be astonishingly close to the numerical results over the 
whole field range. At low field the Doppler shift method 
gives a density of states proportional to £/&a for the d- 
wave case, which corresponds to Volovik's \^B av law. As 
seen already in the single vortex case above, the size of 
the slope turns out to be considerably smaller than in 
the numerical solution, however. For the s-wave case the 
Doppler shift method and the BPT method both yield a 
linear B av field dependence at low fields, while the numer- 
ical solution and the AA method give a y/B^r behavior 
as well. We want to caution, however, that in this field 
range the Abrikosov vortex lattice Eq. (0), which we 
were using here, is not an appropriate groundstate wave- 
function anymore and important corrections from higher 
Landau levels are expect edEJ'C-3 In this field range a fully 
selfconsistea|t,calculation of the vortex lattice becomes 
necessary,E3E3 which will lead to further corrections, like 
for example the Kramer-Pesch effected and the vortex 
core shrinking effect .CJ 
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V. CONCLUSIONS 

We made a detailed comparison of a numerical solution 
of the quasiclassical Eilenberger equations for s- and d- 
wave superconductors in the vortex state with several 
approximate methods. We studied both the single vor- 
tex case and the vortex lattice within a magnetic field 
directed along the c-axis. For the single vortex we found 
that the Doppler shift method for a d-wave gap not only 
fails in the vortex core region, but also misses important 
contributions from core states extending into the nodal 
directions far away from the vortex core. These contribu- 
tions give important corrections to the density of states 
especially at low energies, and thus affect the thermo- 
dynamics of the system. In particular, corrections to 
Volovik's law are found, which are expected to be very 
sensitive to impurity scattering and should be observable 
via systematic impurity studies of the field dependence of 
the specific heat. We expect quantum mechanical effects 
like the Aharanov-Bohm effect to become visible only in 
very clean samples. 

In the vortex lattice there are other approximate 
methods, which are preferred over the Doppler shift 
method. Here, we studied two methods which are valid 



near the upper critical field B C 2'. the method due to 
Brandt, Pesch, and Tewordt and an approximate an- 
alytical method, generalizing a method due to Pesch. 
We showed how the method due to Brandt, Pesch, and 
Tewordt can be derived from the Eilenberger equations 
using the Riccati equations. At low fields both of these 
methods are not able to get the Tomasch resonances, 
which are present in the numerical solution of the Eilen- 
berger equations. However, especially the approximate 
analytical method is impressively close to the numerical 
results at low energies over the whole field range. Since 
this method is also more convenient than the method 
due to Brandt, Pesch, and Tewordt, avoiding a solution 
of an implicit equation and giving analytical results in 
the clean limit, we recommend the use of this method in 
the high field range. 
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